%% initialize
clear all;
close all;
clc

cd  E:\ReplicateBuild\02_code\03_matlab


% diary('mixed_logit')
% diary off
rng(328120)


%%  bring in actual allocations, and then plot against each other (status quo)
T = readtable('E:/03_output/data_fit.csv');
%T=readtable('E:/03_output/data_fit_new_teachers.csv');
fracdis_data=T.p_m2_ma;
black_data=T.black;
aa_data=T.AA;
exp7p_data=T.exp7p;
ca_data=T.CA;


 disp('principal, restricted choice sets')
 load  '../../01_data/05_temp/model/alloc_teach_est_princ_est_restr_baseline'
 
 fracdis=squeeze(jass_iprop_char(:,1,:));
 fracdis_long=fracdis(:);
 
 % absolute advantage of teachers
 AA_iprop=squeeze(jass_iprop_char(:,2,:));
 AA_iprop_long=AA_iprop(:);
 [fracdis_mean,AA_iprop_mean] = bin_data(fracdis_long,AA_iprop_long,20);
 model_fit=fitglm(fracdis_long,AA_iprop_long);
 
 [fracdis_mean_data,AA_mean_data]=bin_data(fracdis_data,aa_data,20);
data_fit=fitglm(fracdis_data,aa_data);

 figure(1)
  plot(fracdis_mean_data, table2array(data_fit.Coefficients(1,1))+table2array(data_fit.Coefficients(2,1))*fracdis_mean_data,'--r')
 hold on
   plot(fracdis_mean, table2array(data_fit.Coefficients(1,1))+table2array(model_fit.Coefficients(2,1))*fracdis_mean,'-k')
%  plot(fracdis_mean, table2array(model_fit.Coefficients(1,1))+table2array(model_fit.Coefficients(2,1))*fracdis_mean,'-k')
 scatter(fracdis_mean_data,AA_mean_data,'r')
 scatter(fracdis_mean,AA_iprop_mean+ table2array(data_fit.Coefficients(1,1))- table2array(model_fit.Coefficients(1,1)),'k','filled')
 hold off
 xlabel('Share disadvantaged students')
 ylabel('Mean absolute advantage of teachers')
 legend('data','model','Location','Southeast')
 saveas(gcf,'../../03_output/figures/AA_validation.png')
 
 
  figure(11)
  plot(fracdis_mean_data, table2array(data_fit.Coefficients(1,1))+table2array(data_fit.Coefficients(2,1))*fracdis_mean_data,'--k')
 hold on
   plot(fracdis_mean, table2array(data_fit.Coefficients(1,1))+table2array(model_fit.Coefficients(2,1))*fracdis_mean,'-k')
%  plot(fracdis_mean, table2array(model_fit.Coefficients(1,1))+table2array(model_fit.Coefficients(2,1))*fracdis_mean,'-k')
 scatter(fracdis_mean_data,AA_mean_data,'k')
 scatter(fracdis_mean,AA_iprop_mean+ table2array(data_fit.Coefficients(1,1))- table2array(model_fit.Coefficients(1,1)),'k','filled')
 text(0.28, 0.0002,'data \rightarrow')
 % text(0.45, 0.0037,'data \rightarrow')
 text(0.665,-0.033,'model \uparrow')
 %annotation('arrow', [0.2 0.22] ,[0.5 0.48]);
 hold off
 xlabel('Share disadvantaged students')
 ylabel('Mean absolute advantage of teachers')
 legend('data','model','Location','Southeast')
 saveas(gcf,'../../03_output/figures/AA_validation_BW.eps')
 
 
  % comparative advantage of teachers
  CA_iprop=squeeze(jass_iprop_char(:,5,:));
 CA_iprop_long=CA_iprop(:);
 [fracdis_mean,CA_iprop_mean] = bin_data(fracdis_long,CA_iprop_long,20);
  model_fit=fitglm(fracdis_long,CA_iprop_long);
 
 [fracdis_mean_data,CA_mean_data]=bin_data(fracdis_data,ca_data,20);
data_fit=fitglm(fracdis_data,ca_data);

 figure(2)
  plot(fracdis_mean_data, table2array(data_fit.Coefficients(1,1))+table2array(data_fit.Coefficients(2,1))*fracdis_mean_data,'--r')
 hold on
  plot(fracdis_mean, table2array(data_fit.Coefficients(1,1))+table2array(model_fit.Coefficients(2,1))*fracdis_mean,'-k')
 scatter(fracdis_mean_data,CA_mean_data,'r')
 scatter(fracdis_mean,CA_iprop_mean+ table2array(data_fit.Coefficients(1,1))- table2array(model_fit.Coefficients(1,1)),'k','filled')
 hold off
 xlabel('Share disadvantaged students')
 ylabel('Mean comparative advantage of teachers')
 legend('data','model','Location','Southeast')
 saveas(gcf,'../../03_output/figures/CA_validation.png')
 
  figure(21)
  plot(fracdis_mean_data, table2array(data_fit.Coefficients(1,1))+table2array(data_fit.Coefficients(2,1))*fracdis_mean_data,'--k')
 hold on
  plot(fracdis_mean, table2array(data_fit.Coefficients(1,1))+table2array(model_fit.Coefficients(2,1))*fracdis_mean,'-k')
 scatter(fracdis_mean_data,CA_mean_data,'k')
 scatter(fracdis_mean,CA_iprop_mean+ table2array(data_fit.Coefficients(1,1))- table2array(model_fit.Coefficients(1,1)),'k','filled')
 text(0.275,0.0034,'data \rightarrow')
 text(0.695,0.0033,'model \downarrow')
 hold off
 xlabel('Share disadvantaged students')
 ylabel('Mean comparative advantage of teachers')
 legend('data','model','Location','Southeast')
 saveas(gcf,'../../03_output/figures/CA_validation_BW.eps')
 
 
     % share teachers with 7+ years experience
 exp7p_iprop=squeeze(jass_iprop_char(:,4,:));
 exp7p_iprop_long=exp7p_iprop(:);
 [fracdis_mean,exp7p_iprop_mean] = bin_data(fracdis_long,exp7p_iprop_long,20);
 
  model_fit=fitglm(fracdis_long,exp7p_iprop_long);
 
 [fracdis_mean_data,exp7p_mean_data]=bin_data(fracdis_data,exp7p_data,20);
data_fit=fitglm(fracdis_data,exp7p_data);

 figure(3)
  plot(fracdis_mean_data, table2array(data_fit.Coefficients(1,1))+table2array(data_fit.Coefficients(2,1))*fracdis_mean_data,'--r')
 hold on
  plot(fracdis_mean, table2array(data_fit.Coefficients(1,1))+table2array(model_fit.Coefficients(2,1))*fracdis_mean,'-k')
 scatter(fracdis_mean_data,exp7p_mean_data,'r')
 scatter(fracdis_mean,exp7p_iprop_mean+ table2array(data_fit.Coefficients(1,1))- table2array(model_fit.Coefficients(1,1)),'k','filled')
 hold off
 xlabel('Share disadvantaged students')
 ylabel('Share teachers with 7+ years of experience')
 legend('data','model','Location','Southeast')
 saveas(gcf,'../../03_output/figures/exp7p_validation.png')
 

 figure(31)
  plot(fracdis_mean_data, table2array(data_fit.Coefficients(1,1))+table2array(data_fit.Coefficients(2,1))*fracdis_mean_data,'--k')
 hold on
  plot(fracdis_mean, table2array(data_fit.Coefficients(1,1))+table2array(model_fit.Coefficients(2,1))*fracdis_mean,'-k')
 scatter(fracdis_mean_data,exp7p_mean_data,'k')
 scatter(fracdis_mean,exp7p_iprop_mean+ table2array(data_fit.Coefficients(1,1))- table2array(model_fit.Coefficients(1,1)),'k','filled')
 text(0.065,0.853,'data \rightarrow')
 text(0.40,0.753,'model \rightarrow')
 hold off
 xlabel('Share disadvantaged students')
 ylabel('Share teachers with 7+ years of experience')
 legend('data','model','Location','Southeast')
 saveas(gcf,'../../03_output/figures/exp7p_validation_BW.eps')
 
 
 % share teachers that are Black
 black_iprop=squeeze(jass_iprop_char(:,3,:));
black_iprop_long=black_iprop(:);
 [fracdis_mean,black_iprop_mean] = bin_data(fracdis_long,black_iprop_long,20);
   model_fit=fitglm(fracdis_long,black_iprop_long);
 
 [fracdis_mean_data,black_mean_data]=bin_data(fracdis_data,black_data,20);
data_fit=fitglm(fracdis_data,black_data);

 figure(4)
  plot(fracdis_mean_data, table2array(data_fit.Coefficients(1,1))+table2array(data_fit.Coefficients(2,1))*fracdis_mean_data,'--r')
 hold on
  plot(fracdis_mean, table2array(data_fit.Coefficients(1,1))+table2array(model_fit.Coefficients(2,1))*fracdis_mean,'-k')
 scatter(fracdis_mean_data,black_mean_data,'r')
 scatter(fracdis_mean,black_iprop_mean+ table2array(data_fit.Coefficients(1,1))- table2array(model_fit.Coefficients(1,1)),'k','filled')
 hold off
 xlabel('Share disadvantaged students')
 ylabel('Share Black teachers')
 legend('data','model','Location','Southeast')
 saveas(gcf,'../../03_output/figures/black_validation.png')
 
  figure(41)
  plot(fracdis_mean_data, table2array(data_fit.Coefficients(1,1))+table2array(data_fit.Coefficients(2,1))*fracdis_mean_data,'--k')
 hold on
  plot(fracdis_mean, table2array(data_fit.Coefficients(1,1))+table2array(model_fit.Coefficients(2,1))*fracdis_mean,'-k')
 scatter(fracdis_mean_data,black_mean_data,'k')
 scatter(fracdis_mean,black_iprop_mean+ table2array(data_fit.Coefficients(1,1))- table2array(model_fit.Coefficients(1,1)),'k','filled')
 text(0.15,0.215,'data \rightarrow')
 text(0.52,0.33,'model \rightarrow')
 hold off
 xlabel('Share disadvantaged students')
 ylabel('Share Black teachers')
 legend('data','model','Location','Southeast')
 saveas(gcf,'../../03_output/figures/black_validation_BW.eps')
 
 
  
%% plot distribution (christina's figure)
clear
% 0. first best

load  '../../01_data/05_temp/model/alloc_teach_fracdis_princ_va_unrestr_baseline'
%load '../../01_data/05_temp/thetahat'

firstbest_school_store=firstbest_school;

clear  stable_teach_propose_store teach_school_same_store 

% 1. status quo:  alloc_teach_princ_restr
load  '../../01_data/05_temp/model/alloc_teach_est_princ_est_restr_baseline'
clear i* jass* teach_school* 

stable_teach_propose_teach_princ_restr=stable_teach_propose;
stable_school_propose_teach_princ_restr=stable_school_propose;

clear stable_school_propose  stable_teach_propose firstbest_school 

% 2. timing:  alloc_teach_est_princ_est_unrestr_

load  '../../01_data/05_temp/model/alloc_teach_est_princ_est_unrestr_baseline'
clear i* jass* teach_school* random* jrank*


stable_teach_propose_teach_princ_unrestr=stable_teach_propose;

clear stable_school_propose  stable_teach_propose  firstbest_school 

% 3. principal max VA, restricted choice sets

load  '../../01_data/05_temp/model/alloc_teach_est_princ_va_restr_baseline'
clear i* jass* teach_school*

stable_teach_propose_teach_est_princ_va_restr=stable_teach_propose;

clear stable_school_propose  stable_teach_propose firstbest_school 

% 4a. teacher max n dis, restricted choice sets 

load  '../../01_data/05_temp/model/alloc_teach_ndis_princ_est_baseline'
clear i* jass* teach_school* 

stable_teach_propose_teach_ndis_princ_restr=stable_teach_propose;

clear stable_school_propose  stable_teach_propose firstbest_school

% 5.  teacher max frac dis, principal VA max, timing restricted

load  '../../01_data/05_temp/model/alloc_teach_ndis_princ_va_baseline'
clear i* jass* teach_school* 

stable_teach_propose_teach_ndis_princ_va_restr=stable_teach_propose;

clear stable_school_propose  stable_teach_propose firstbest_school 


% figure for the paper
figure(5)
scatter(mean(stable_teach_propose_teach_princ_restr(:,3)), mean(stable_teach_propose_teach_princ_restr(:,4)),100,'ko')
set(gca,'FontSize',15)
hold on
scatter(mean(stable_school_propose_teach_princ_restr(:,3)), mean(stable_school_propose_teach_princ_restr(:,4)),100,'ko')
scatter(mean(stable_teach_propose_teach_princ_unrestr(:,3)), mean(stable_teach_propose_teach_princ_unrestr(:,4)),100,'ko')
scatter(mean(stable_teach_propose_teach_est_princ_va_restr(:,3)), mean(stable_teach_propose_teach_est_princ_va_restr(:,4)),100,'ko')
scatter(mean(stable_teach_propose_teach_ndis_princ_restr(:,3)), mean(stable_teach_propose_teach_ndis_princ_restr(:,4)),100,'ko')
scatter(mean(stable_teach_propose_teach_ndis_princ_va_restr(:,3)), mean(stable_teach_propose_teach_ndis_princ_va_restr(:,4)),100,'ko')
scatter( mean(firstbest_school_store(:,3)), mean(firstbest_school_store(:,4)),100,'ko')
plot( (-0.14:0.01:0.04)',(-0.14:0.01:0.04)','--k','Linewidth',1)
hold off
ylim([-0.14 0.04])
xlim([-0.14 0.05])
xlabel('Advantaged students (VA)')
ylabel('Disadvantaged students (VA)')
saveas(gcf,'../../03_output/figures/dist_everything.png')

figure(51)
scatter(mean(stable_teach_propose_teach_princ_restr(:,3)), mean(stable_teach_propose_teach_princ_restr(:,4)),100,'ko')
set(gca,'FontSize',15)
hold on
scatter(mean(stable_school_propose_teach_princ_restr(:,3)), mean(stable_school_propose_teach_princ_restr(:,4)),100,'ko')
scatter(mean(stable_teach_propose_teach_princ_unrestr(:,3)), mean(stable_teach_propose_teach_princ_unrestr(:,4)),100,'ko')
scatter(mean(stable_teach_propose_teach_est_princ_va_restr(:,3)), mean(stable_teach_propose_teach_est_princ_va_restr(:,4)),100,'ko')
scatter(mean(stable_teach_propose_teach_ndis_princ_restr(:,3)), mean(stable_teach_propose_teach_ndis_princ_restr(:,4)),100,'ko')
scatter(mean(stable_teach_propose_teach_ndis_princ_va_restr(:,3)), mean(stable_teach_propose_teach_ndis_princ_va_restr(:,4)),100,'ko')
scatter( mean(firstbest_school_store(:,3)), mean(firstbest_school_store(:,4)),100,'ko')
plot( (-0.14:0.01:0.04)',(-0.14:0.01:0.04)','--k','Linewidth',1)
hold off
ylim([-0.14 0.04])
xlim([-0.14 0.05])
xlabel('Advantaged students (VA)')
ylabel('Disadvantaged students (VA)')
saveas(gcf,'../../03_output/figures/dist_everything_BW.eps')

